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SHOCK  AND  DAMAGE  IN  MARINE  COMPOSITE  STRUCTURES 


1.0  INTRODUCTION 

The  use  of  composite  structures  for  present  and  future  naval  applications  represents  an 
important  development.  Composites  offer  significant  advantages  over  the  traditional  metals  due 
to  their  strength,  stiffness  and  lightweight  characteristics.  However,  the  behavior  of  these 
composite  structures,  particularly  under  highly  transient  shock  loadings,  is  not  well  understood  at 
the  present  time.  Material  and  structural  failure  models,  subjected  to  dynamic  loading 
environments  are  not  currently  well  developed  especially  for  thick,  polymer  matrix,  fiber 
reinforced,  composite  materials.  Also,  the  experimental  data  base  for  these  materials  in  naval 
type  structures  is  very  sparse.  References  [1-3]  are  indicative  of  some  of  the  experimental  work 
that  has  been  performed  on  composite  structures  subjected  to  underwater  shock  environments. 

Composite  structures  represent  a  radical  departure  from  the  traditional  ductile, 
homogeneous,  isotropic  metal  structures.  Under  severe  loadings,  metals  deform  plasticity  due  to 
slip  along  shear  lines.  Current  analysis  methods,  which  are  rooted  in  continuum  mechanics,  have 
been  developed  and  applied  based  upon  traditional  metal  type  of  structures  where  failures  usually 
have  been  more  global  in  nature  due  to  the  metal  ductility.  The  use  of  high  strength  metals  in 
naval  structures  has  lead  to  more  brittle  types  of  problems,  such  as  weld  cracks.  But  modem 
finite  element  techniques  and  plasticity  constitutive  models  can  still  be  widely  applied  due  to  the 
homogeneous  and  nearly  isotropic  nature  of  the  metals. 

Fiber  reinforced  polymer  composites  are  very  heterogeneous  materials  that  combine  high 
performance  fibers  in  a  viscoelastic  matrix.  There  are  several  levels  of  heterogeneities  to 
consider:  between  layers,  between  the  fiber  and  matrix,  and  also  heterogeneities  in  the  form  of 
voids  arising  from  processing  the  composite.  Thus  composites  are  more  difficult  to  characterize 
and  model  with  the  standard  “deterministic”,  continuum  mechanics  based,  finite  element 
programs  currently  available. 

In  addition,  composite  structures  tend  to  be  highly  dispersive  to  propagating  waves,  and 
deform  nonlinearly  under  severe  load  due  to  the  development  of  networks  of  micro-cracks. 
Because  of  all  these  features,  the  shock  response  and  the  development  of  damage  in  composites 
is  not  a  well  understood  phenomena. 
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Continuum  damage  mechanics,  see  for  instance  [4-14],  is  in  general  a  developing  field  of 
research  that  is  not  yet  mature.  Damage  theories  have  progressed  further  for  metal  structures  due 
to  their  more  homogeneous,  isotropic  and  ductile  nature.  The  abilitiy  to  model  and  predict 
localization  and  shear  banding  as  well  as  the  propagation  of  large  cracks  has  certainly  progressed 
in  recent  years.  Plasticity  assumptions  which  allow  the  use  of  the  scalar  variables,  such  as 
effective  stress  and  strain,  are  responsible  for  much  of  this  success  in  metal  structures. 

Much  of  the  modem  concept  of  damage  mechanics  is  based  upon  the  idea  of  area  reduction. 
Kachanov  [15]  in  1958  first  purposed  a  phenomenological  theory  of  creep  damage  by  introducing 
a  scalar  damage  variable  D  to  characterize  the  material  degradation.  Later  Rabotnov  [16] 
interpreted  D  as  representing  the  fraction  of  damaged  material.  The  effective  resisting  area  of 
the  local  material  degrades  initially  from  A  to  A  due  to  damage  where: 

A  =  (1  -  D)A  (1.1) 

In  one  dimension  (ID),  the  effective  stress  is  defined  as: 

a  =  o/(\-D )  (1.2) 

This  concept  of  local  material  damage  and  effective  resisting  area  has  been  extended  to  include 
vector  and  higher  order  tensor  definitions  of  damage  more  applicable  to  orthotropic  and 
anisotropic  material.  The  literature  in  continuum  damage  is  quite  extensive,  see  for  instance  [4- 
11]  and  [24-26]  for  more  of  a  discussion  regarding  metals  and  plasticity,  and  [12-14]  for 
discussion  of  brittle  and  composite  materials. 

Although  such  concepts  as  local  damage  described  by  D ,  which  may  be  represented  as  a 
scaler,  vector  or  higher  order  tensor,  are  quite  useful  one  must  always  remember  that  composites 
are  a  very  heterogeneous  media  in  general.  Continuum  based  models  are  therefore  particularily 
difficult  to  apply  given  the  statistical  nature  of  the  media.  The  linking  between  scales,  micro- 
meso-macro  (structural)  is  also  extremely  important  but  at  present  not  well  understood.  Issues 
associated  with  homogenization  (see  for  instance  [17])  must  be  adequately  addressed.  Because 
of  these  factors,  much  of  what  has  been  and  is  being  developed  in  the  continuum  damage  field  is 
limited  to  specific  materials  and  structures  and  relatively  narrow  in  application. 
Phenomenological  approaches  have  been  prevalent. 

One  very  important  exception  to  this  restriction  is  the  dissipated  energy  approach  developed  by 
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Mast  et  al  in  [18].  In  this  work,  the  authors  couple  the  calculation  of  the  damage  in  a  test 
specimen  directly  to  the  computation  of  the  dissipated  energy.  An  assumption  is  made  that  the 
material  behaves  as  indicated  in  Figure  1.1,  where  where  the  stress  strain  curve  follows  path 
ABC  for  loading  and  by  path  CA  for  unloading.  The  area  between  curves  ABC  and  CA 
represents  the  dissipated  energy  of  the  material.  Using  the  experimental  data,  a  dissipated  energy 
density  function  <p  is  obtained  for  the  composite  material  via  interpolation.  The  stress  c  is 

postulated  to  be  a  function  of  strain  and  a  single  state  variable  £  or: 

g  =  C(^,e)  (1.3) 

where  C  is  the  constitutive  tensor.  Next  the  work  potential: 

^(e)  =  3>(e)  +  <p(e)  (1.4) 

where  d>(£)  is  equal  to  the  energy  density  recovered  if  the  material  were  to  unload  linearly  from 
its  present  state.  The  constitutive  tensor  is  then  determined  from: 

5  =  ¥*¥(£,£)  (1.5) 

The  analytical  and  experimental  approach  in  [18]  is  limited  to  static  or  quasi-static  conditions, 
can  not  detect  through  thickness  delaminations,  and  defines  damage  only  in  a  scalar  sense.  Thus, 
directional  damage  can  not  be  differentiated.  However,  it  has  been  highly  successful  and  is  a 
very  important  contribution  to  the  field  of  damage  mechanics  because  it  directly  couples  a  series 
of  experiments  and  loading  conditions  for  composite  materials  to  the  calculation  of  the 
constitutive  tensor  and  damage.  This  work  has  the  added  advantage  that  unlike  other  approaches, 
it  does  not  require  the  estimation  of  several  additional  parameters  that  are  good  only  for  a  narrow 
range  of  material  properties. 

In  general  in  this  report,  we  shall  emphasize  the  evolving  damage  approach  and  results 
developed  in  [19-22],  In  section  2,  ID  and  2D  continuum  damage  models  are  reviewed  and 
compared  to  the  experimental  results  [23]  obtained  from  impact  tests  on  thick  composite  circular 
plates.  In  section  3,  ID  parametric  studies  of  graphite/peek  plates  subjected  to  underwater  shock 
(UNDEX)  are  discussed.  Approximate  ranges  for  incipient  damage  and  extensive  damage  are 
indicated  for  the  weight  of  the  charge  and  its  standout  distance. 

In  section  4,  the  3D  continuum  damage  theory  developed  in  [22]  is  summarized.  A  preliminary 
version  of  the  3D  damage  model  has  been  programmed  into  the  commercial  finite  element 
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program  ABAQUS/EXPLICIT  but  not  yet  tested.  In  general  as  discussed  in  section  4  and  [22], 
the  modeling  of  evolving  3D  damage  in  composite  plates  and  shells  represents  an  extremely 
difficult  task  to  do  efficiently  and  accurately.  Several  hurdles  remain  to  overcome  before  robust 
3D  capabilities  will  be  available  for  general  use. 

Finally  in  section  5,  a  summary  and  some  conclusions  are  given,  from  the  perspective  of  work 
done  here  at  NRL,  regarding  the  state  of  the  art  of  modeling  the  shock  response  and  damage  in 
composites. 

1.1  Through  thickness  modeling  and  damage 

As  we  have  seen,  composite  structures  differ  significantly  from  metal  structures  in  their 
response  to  loading  environments.  Composite  plates  and  shells  are  typically  thicker  than 
comparable  metal  structures.  Couple  this  with  the  lower  compliances  of  their  matrix  components, 
and  fiber  reinforced  composites  typically  display  much  more  transverse  shear  effects  in  general. 
Current  finite  element  codes  can  account  only  in  a  crude  fashion  for  transverse  shear  effects. 

Fiber  reinforced  composites  are  also  more  heterogeneous  and  anisotropic  than  typical  structural 
metals.  The  layering  of  various  ply  combinations  produces  a  complex  arrangement  that  can  only 
be  approximate  in  an  average  sense.  Without  some  sort  of  homogenization  procedure,  one  is 
forced  to  attempt  a  detailed  modeling  throught  the  thickness  direction  of  the  plate  or  shell 
structure  using  solid  elements.  This  is  typically  completely  unfeasible  given  the  thickness  to 
length  ratios  of  typical  plate/shell  structures,  which  would  require  large  number  of  solid  elements. 

Through  thickness  (or  delamination)  damage,  especially  in  dynamic  loading  environments,  is 
much  more  of  a  concern  in  composite  structures  than  it  is  in  metal  structures.  This  type  of 
damage  is  also  very  difficult  to  accurately  model.  The  standard  displacement  finite  element  codes 
that  are  popular  today  do  not  contain  plate  or  shell  elements  which  can  account  for  through 
thickness  (normal  direction)  deformation  and  stress  even  in  a  crude  fashion.  These  types  of 
plate/shell  elements  simply  do  not  exist. 

To  account  for  through  thickness  effects  using  the  currently  available  commercial  finite  element 
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method  (FEM)  technology  one  has  no  other  choice  but  to  use  solid  elements  -  thus  producing 
very  large  models  to  solve.  In  sections  2,  3  and  4,  more  will  be  said  about  modeling  through  the 
thickness. 

1.2  Dynamic  Versus  Static  Behavior 

Composite  structures,  especially  laminated  plates  and  shells,  comprise  a  very  heterogeneous 
media,  whose  dynamic  response  can  be  exceedingly  complex.  If  the  wavelength  of  the  loading 
and  the  response  of  the  material  is  very  long  compared  to  the  scale  of  the  inhomogeneity,  then  the 
material  response  is  governed  by  effective  properties  of  the  equivalent  homogenized  media. 
However,  if  the  composite  structure  is  subjected  to  shock  environments,  then  the  wavelengths  of 
the  loading  and  response  of  the  structure  will  be  much  shorter.  In  this  case,  the  characteristic 
dimensions  of  the  heterogeneous  media  become  much  more  important.  The  interfaces  between 
the  material  phase  cause  wave  reflection  and  refraction.  The  energy  is  thus  spread  or  “dispersed’ 
over  many  wavelengths. 

Dispersion  can  also  be  present  in  composite  materials  due  to  the  viscoelastic  nature  of  the 
matrix.  Viscoelasticity  represents  a  history  dependent  nonlinear  response  that  is  difficult  to 
include  even  approximately  in  numerical  models.  The  response  of  viscoelastic  material  is  also 
usually  frequency  dependent.  Couple  this  to  the  general  anisotropic  and  heterogeneous  nature  of 
a  laminated  structure,  and  the  inclusion  of  viscoelasticity  effects  into  numerical  models  becomes 
even  more  difficult  to  accurately  account  for. 

In  [21],  Nemes  and  Randles  discuss  the  inclusion  of  viscoelasticity  into  a  ID  damage  model. 
They  obtain  an  excellent  fit  to  experimental  data  [23]  for  the  impact  of  a  flyer  plates  on 
axisymmetric  graphite/peek  composite  plates.  Again  this  work  is  ID,  in  the  through  thickness 
direction. 

Viscoelasticity  and  material  damping  remain  a  significant  problem  to  overcome  for  numerical 
models  of  even  linear  composite  structures.  Current  finite  element  modeling  approaches  appear 
to  use  various  forms  of  structural  damping  such  as  Rayleigh  to  try  to  account  for  the  damping 
caused  by  the  viscoelasticity.  For  models  that  include  evolving  damage,  other  than  [21,22] 
explicit  viscoelasticity  effects  have  been  largely  ignored. 
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With  respect  to  metal  structures,  Nemes  et  al.  in  a  series  of  papers  [24-26]  developed  a 
modified  form  of  the  viscoplastic  damage  constitutive  theory  of  Perzyna  [27]  and  applied  it  to 
model  high  strain  rate  behavior  occuring  in  plate  impact  spall  fracture.  Although  they  were  quite 
successful  and  obtained  excellent  correlation  to  experimental  data,  this  type  of  approach,  which  is 
phenomenological  based,  requires  extensive  experimental  data  to  even  estimate  the  several 
material  and  rate  parameters  that  are  needed.  The  range  of  applicability  of  this  viscoplastic 
damage  approach  is  somewhat  narrow,  and  of  course  restricted  to  metals. 

In  a  dynamic  loading  situation,  damage  to  the  structure  evolves  in  a  very  short  period  of  time  . 
Many  of  the  thermodynamically  based  approaches  are  simply  not  completely  valid  because  the 
structure  is  not  in  equilibrium.  Damage  occurs  very  quickly. 

As  mentioned  in  section  1.1,  through  thickness  or  delamination  damage  for  a  composite 
structure  is  an  important  consideration  especially  for  a  dynamic  load  such  as  impact  or  UNDEX. 
Composites  are  much  more  susceptible  to  this  type  of  damage  because  they  are  very 
heterogeneous,  are  ususally  thicker,  have  a  fairly  low  modulus  for  the  matrix,  and  a  lower  wave 
speed  (especially  throught  the  thickness)  than  comparable  metal  structures.  For  instance  under  an 
impact  loading  to  a  plate  or  shell  structure,  the  compression  wave  passes  through  the  thickness 
and  is  reflected  off  the  backface  as  a  tension  wave.  If  the  unloading  or  release  wave,  which  is  also 
tensile,  meets  this  reflected  wave  delamination  can  occur  in  a  composite  if  the  tensile  strength  of 
the  matrix  is  exceeded.  Because  the  composite  is  typically  thicker  and  has  a  lower  wave  speed 
than  its  metal  counterpart,  the  release  wave  has  more  time  and  distance  to  meet  the  reflected  wave 
and  set  up  this  situation  for  the  composite. 
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Figure  1 . 1  Assumed  material  behavior 
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2.  ID  AND  2D  CONTINUUM  DAMAGE 


In  this  section,  we  first  summarize  the  NRL  2D  continuum  damage  model  (CDM)  developed  in 
[19,20],  This  is  followed  by  a  comparison  of  ID  and  2D  models  to  available  experimental  data 
[23]  in  subsections  2.2  and  2.3. 

2. 1  Summary  of  the  NRL  2D  Continuum  Damage  Model 

References  [19,20]  develop  a  2D  continuum  damage  model  (CDM)  that  employs  an  evolving 
vector  damage  theory  for  transversely  isotropic  materials.  This  theory  is  intended  for  application 
to  highly  transient  loading  situations  such  as  impact  and  underwater  shock.  Figure  2.1  indicates  a 
typical  laminated  composite.  An  assumption  of  transverse  isotropy  is  made.  The  matrix  damage 
is  characterized  by  a  vector  or  [14]: 

V  =  V,el+V!e!+V!eJ  (2.1) 

where  e,  i=l,2,3  are  the  unit  normals  along  the  material  coordinate  directions  shown  in  Figure 
2  .1.  In  this  phenomenological  treatment,  the  damage  magnitude  is  interpreted  simply  in  terms  of  a 
fractional  reduction  in  certain  of  the  elastic  properties.  The  V,  component  represents  a  network 
of  matrix  cracks  in  the  2-3  plane  which  induce  delaminations  throught  the  thickness  direction  of 
the  laminate.  The  V2  and  V3  components  of  damage  are  caused  by  matrix  cracks  which  may 

traverse  or  lie  between  the  reinforcing  fibers,  but  do  not  break  the  fibers.  Since  a  transverse 
isotropic  assumption  is  made,  the  inplane  components  of  the  damage  vector  V2and  V3  are 
combined  into  the  scalar  Vs  by: 

VS=(V2+V2),/2  (2.2) 

Thus,  the  inplane  damage  Vsis  characterized  by  a  scalar  with  a  general  direction  that  is 
perpendicular  to  the  V,  damage. 

The  evolving  damage  components  V,  and  Vs  vary  from  0  (undamaged)  to  1  (complete  failure 

at  the  material  point).  The  damage  components  are  assumed  by  Nemes  and  Randles  [19,20]  to 
reduce  the  elastic  properties  of  the  composite  at  the  material  point  in  the  transversely  isotropic 
media  as  follows: 
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En  =(1- V,2)E°, 

E22=(l-a,Vs2)E°2 

Gl2  =  (1  -  V,2)(l  -  a2Vs2)Gf2  (2.3) 

v12=(l-V2)(l-a3V2K2 
v23=(l-a4V2)v°3 

where  E,  G,  and  v  denote  the  directional  elastic  moduli,  shear  moduli  and  Poisson’s  ratios.  The 
superscript  “0”  denotes  virgin  properties  and  the  fractions  CXa^l,  i=l  -  4  are  included  to 
prevent  a  complete  loss  of  material  integrity  as  Vs  — >  1. 

Damage  evolution  and  rate  dependency  is  introduced  in  the  2D  Nemes  and  Randles  CDM 
approach  by  assuming  a  dependency  on  the  current  state  of  damage,  some  stress  above  a  current 
threshold,  and  material  properties  controlling  damage  evolution  rates.  Both  V,  and  Vs  types  of 

damage  are  assumed  to  be  governed  by  a  scalar  threshold  function  F  (analogous  to  a  yield 
function  in  plasticity)  of  the  form: 

F(c,  f(V))  <  0  for  no  growth  and  (2.4) 

>  0  for  damage  growth 

where  o  is  the  current  stress  tensor,  f  is  an  array  of  current  threshold  parameters  which  is  a 
function  of  V ,  the  current  damage  vector. 

The  threshold  function  F  is  assumed  by  Nemes  and  Randles  to  be  similar  in  form  to  the  Mohr- 
Coulomb  [28]  type  of  yield  criteria  used  in  plasticity.  This  allows  the  strength  of  the  material  to 
be  dependent  on  the  tension.  The  form  of  F  is  assumed  to  be: 

F(o,f)  =  (l  +  (T/f3)1/2-(f1-o)/f2  (2.5) 

In  (2.5),  x  and  o  represent  effective  types  of  shear  and  tension  stresses  and  differ  in  form  for  V, 
or  Vs  type  of  damage  evolution.  The  parameters  f,  1=1, 2, 3  are  related  to  specific  growth 
threshold  strengths  (oG  and  xG)  and  the  Coulomb  friction  tangent  (pGl2  as: 

°gi  =  fi  -  f2  (tension  threshold) 

tgi  =  I3  •  ((^  /  f2)2  - 1)1/2  (shear  threshold)  (2.6) 

^gi  =  I3 1  f2 
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For  the  V,  component  of  the  damage,  a  and  x  in  (2.5)  are  assumed  to  be  of  the  form: 

a  =  on  (tension  only) 

x  =  o,j  (2.7) 

The  growth  threshold  strengths  and  Coulomb  friction  tangent  in  (2.6)  are  then  postulated  to 
depend  on  the  V,  damage  as: 

^gi  =  (1  —  V,  )oG10 

*oi=<l-V,:)tolo  (2.8) 

<Pg«  =  ^GIO  +  Vj  (<PgU  -  ^Gio) 

where  the  “0”  subscript  denotes  virgin  (V,  =0)  threshold  properties.  Equations  (2.8)  are 
substituted  into  (2.6)  and  solved  for  f, ,  f2  and  f3 .  A  similar  approach  is  used  to  determine  f,  for 
i=l,2,3  for  Vs  damage  -  see  [19], 

The  evolution  for  V,  and  Vs  damage  completes  the  CDM  constitutive  description.  The  rate 
forms  for  V,  and  Vs  damage  used  are: 

^  =  (d,  /Ooi,)1'  /(n,(l -  V,!)  (2.9a) 

at 

^■  =  (d,/o0M)-/n,  (2.9b) 

at 

where  n, ,  r|,  and  ns,  qs  are  empircally  based  material  constants,  while  d,  and  ds  represent 
distances  from  the  threshold  (yield)  surfaces  for  V,  and  Vs  damage.  Thus,  if  d,  or  d2  is  zero, 
the  rate  of  the  damage  growth  will  also  be  zero.  Note  in  (2.9a)  as  V,  ->  1 ,  an  acceleration  to  a 
castastrophic  failure  occurs,  while  in  (2.9b)  no  explicit  dependence  on  Vs  is  assumed  for 
dVs  /  dt . 

Overall,  the  2D  CDM  approach  developed  by  Nemes  and  Randles  [19-20]  is  very  appealing  and 
has  had  some  success.  However,  like  all  phenomenological  theories,  it  requires  several  material 
parameters  that  can  only  be  estimated  through  pertinent  experiments  on  specific  laminated 
composites.  In  addition,  this  appproach  like  so  many  others  does  not  adequately  address  the 
explicit  viscoelastic  nature  of  the  matrix  in  the  plane  of  the  fibers.  In  [21,22],  viscoelasticity  is 
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added  to  a  ID  version  of  the  CDM,  and  good  agreement  is  obtained  with  the  experimental  data 
for  the  impacted  composite  plates  [23],  Sections  2.2  and  2.3  contain  further  discussion. 
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2.2  COMPARISONS  OF  ID  MODELS  TO  EXPERIMENTAL  RESULTS 

Impact  experiments  were  conducted  on  axisymmetric  graphite/peek  and  graphite/epoxy  plates 
for  the  U.S.  Naval  Research  Lab  on  the  Phillips  Lab  Impact  Facility  102  mm  gas  gun,  using 
standard  plate-impact  techniques.  Two  types  of  flyer-plate  impact  experiments  were  conducted  - 
see  [23],  The  first  were  transmitted  wave  tests,  while  the  second  set  of  experiments  were 
spallation  tests.  The  flyer  consisted  of  axisymmetric  polymethylmethacrylate  (PMMA)  plates. 

One  dimensional  (ID)  numerical  simulations,  using  the  WONDY  [29]  finite  difference  computer 
program,  for  these  sets  of  experiments  are  reported  in  [21]  and  [22],  the  inclusion  of 
viscoelasticity  effects  in  the  through  thickness  direction  produced  ID  numerical  results  that 
compared  very  well  to  the  experimental  test  data.  In  these  detailed  computer  models,  hundreds  of 
difference  points  or  zones  were  used  to  model  through  the  thickness  of  the  flyer  and  plate  in  order 
to  accurately  describe  the  shock  structure.  Figure  2.2.1  is  taken  from  [21]  and  it  compares  the 
computed  and  measured  particle  velocities  of  the  backface  of  the  plates  in  the  spall  experiments 
(referred  to  as  shots).  In  general,  the  comparison  is  very  good.  In  Figure  2.2.2,  again  from  [21], 
computed  V,  damage  profiles  through  the  plates  are  indicated  for  various  experiments  (shots).  A 
value  near  1.0  indicates  total  failure  at  that  material  point.  These  predictions  compare  reasonably 
well  to  observed  damage  (see  [23])  and  ultrasonic  examination  [21], 

ID  numerical  simulations  with  the  NRL  2D  CDM  have  also  been  investigated  using  the  finite 
element  codes  PR0NT02D  and  more  recently  ABAQUS/EXPLICIT  -  see  [22],  In  these  2D 
models  of  ID  behavior,  explicit  viscoelasticity  has  not  been  included  through  the  thickness,  and 
the  models  are  much  coarser  than  the  very  refined  WONDY  models.  Even  so,  these  results 
compare  quite  well  in  general  to  the  WONDY  runs  and  experimental  data. 
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Figure  2.2.1  Comparison  of  computed  and  measured  particle  velocity  in  spall  experiments 

( . computed, _ experiment):  (a)  shot  3360;  (b)  shot  3357;  (c)  shot 

3362;  (d)  shot  3358;  (e)  shot  3368. 
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2.3  COMPARISONS  OF  2D  MODELS  TO  EXPERIMENTAL  IMPACT  RESULTS 


Coarse  2D  finite  element  models,  using  the  PR0NT02D  program,  were  constructed  in  [22]  and 
compared  well  to  the  experimental  impact  tests  in  [23],  As  indicated  in  section  2.2,  the  2D  NRL 
continuum  damage  model  (described  in  section  2.1)  has  also  been  implemented  into  the  explicit 
finite  element  code  ABAQUS/EXPLICIT,  and  currently  is  available  on  the  NRL  Cray-EL. 

Before  discussing  specific  models,  some  comments  concerning  the  modeling  of  2D  and  of 
course  3D  structures  are  in  order.  In  our  ID  simulations,  hundreds  of  finite  diffemce  points 
(WONDY  analysis)  or  finite  elements  (PR0NT02D  analysis)  were  employed  through  the  plate 
thickness  to  accurately  capture  the  propagation  of  the  shock  wave  and  the  insuing  damage.  For 
2D  analysis  this  would  translate  into  hundreds  of  elements  through  the  thickness.  Due  to  aspect 
ratio  considerations,  this  level  of  refinement  would  require  tens  of  thousands  of  elements  to  model 
just  a  2D  plate  and  flyer  -  which  is  completely  impractical.  In  addition,  2D  analyses  require  much 
longer  times  than  ID  runs  due  to  flexural  and  shear  waves  that  develop  in  the  plates. 

Thus  in  our  limited  number  of  2D  models  we  have  analyzed  here  at  NRL,  our  purpose  has  not 
been  to  accurately  capture  the  propagation  and  structure  of  the  shock  wave,  along  with  the 
flexural  and  shear  waves.  But  rather  we  have  limited  ourselves,  due  primarily  to  computer  costs, 
to  coarser  models  in  which  to  predict  overall  damage  patterns  in  the  plates  and  estimates  of  peak 
backface  velocities. 

Figure  2.3.1  indicates  a  2000  element  axisymmetricmodel  of  a  PMMA  flyer  and  a  graphite/peek 
plate  using  ABAQUS/EXPLICIT.  A  30x50  grid  of  2D,  4  node  elements  were  used  for  the  plate  - 
30  elements  through  the  thickness.  A  grid  of  10x50  elements  was  employed  to  discretize  the 
flyer.  In  the  2D  simulations,  which  were  run  on  the  NRL  Cray-EL,  the  flyer  was  given  an  initial 
velocity  of  93  m/sec,  which  corresponds  to  shot  3357  in  the  ID  runs  in  Figures  2.2.1  and  2.2.2, 
and  the  experimental  test  [23]. 

In  Figure  2.3.2,  the  predicted  backface  velocity  of  the  plate  is  indicated.  Although  the  shape  of 
this  curves  does  not  compare  that  well  to  Figure  2.2.1b,  this  coarse  2D  model  does  predict  a 
maximum  particle  velocity  of  approximately  77  m/sec,  which  is  very  close  to  both  the  refined 
WONDY  prediction  and  the  measured  velocity.  Figure  2.3.3  described  the  predicted  through 
thickness  VI  damage  along  the  centerline  of  the  axisymmetric  graphite/peek  plate  for  the  2000 
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element  ABAQU S/EXPLICIT  model.  This  result  compares  reasonably  well  to  Figure  2.2.2b  and 
the  refined  WONDY  model. 

In  general,  we  are  encouraged  by  these  results  from  coarse  2D  models  with  the  NRL  2D 
continuum  damage  theory.  Applications  to  other  composite  plates  as  well  as  parametric  studies 
can  be  done  with  this  capability,  which  has  been  programmed  as  a  user  written  subroutine  for 
ABAQUS/EXPLICIT. 
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Figure  2.3.1  2000  element  axisymmetric  model  of  a  PMMA  flyer  impacting  a 
graphite/epoxy  plate. 


backface  (centerline)  velocity  (m/sec) 


Delamination  (VI)  damage 


Relative  distance  through  the  thickness  of  the  plate  (0  is  the  top  and  1  is  the  backface) 


Figure  2.3.3  Delamination  (VI )  damage  through  the  centerline  of  the  plate 
for  the  2000  element  model. 
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3.  FLUID-STRUCTURE  INTERACTION  MODELS 


In  this  section  ID  results  from  the  modeling  of  graphite/peek  plates  subjected  to  underwater 
shock  are  examined.  The  evolving  damage  model  developed  by  Nemes  and  Randles  [19,20]  has 
been  employed. 

3.1  ID  Models 

In  [22],  a  ID  version  of  the  2D  evolving  damage  model  of  Nemes  and  Randles  was  extended  to 
include  dispersion  effects.  A  series  of  simulations  were  performed  with  the  WONDY  [29]  finite 
difference  computer  program  to  predict  the  response  of  composite  graphite/peek  plates  from  2  to 
8  inches  thick.  The  plates  were  subjected  to  simulated  near  field  underwater  explosions.  The 
intent  of  these  analyses,  which  are  documented  in  [22],  was  to  assess  the  extent  of  spallation 
damage  for  various  charge  weights  and  ranges,  and  to  determine  the  approximate  boundaries  of 
incipient  and  complete  spallation.  The  history  of  a  TNT  explosive  pulse  passing  through  a  fixed 
spatial  point  can  be  expressed  as  [30]: 

P  =  Poe~‘/e  (3.1) 

where  p0  is  the  initial  pressure  and  0  the  exponential  decay  constant.  Using  Aaron’s  law  p0  and 
0  can  be  expressed  as: 

W1/3 

p0  =  0.519(—  .)*■'>.  (3.2) 

0  =  0.925x10^ 

in  which  p0  is  given  in  kilobars,  0  in  seconds,  R  is  the  range  in  meters  (m),  and  W  is  the  charge 

mass  in  kilograms  (kg).  The  value  of  W  in  (3.2)  for  this  study  ranged  from  5  kg  to  30  kg,  with  R 
varying  from  .5  to  1 .5  m. 

For  seawater,  the  following  material  properties  (SI  units)  were  employed  in  the  modeling: 

pw  =  1000  kg/m3  (3.3) 

cw  =  1500  m/s 
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sw  =  1.75 

where  pw  is  the  density  and  cw  is  the  water  speed.  The  variable  sw  comes  from  a  linear  shock 
velocity  us  to  particle  velocity  up  relationship  of  the  form: 

us=cw+sw-up  (3.4) 

The  material  constants  used  for  the  graphite/peek  plates  were: 
bulk  properties 

p=  1579  kg/m3  (3.5a) 

c0  =  3000  m  /  s  (wave  speed  through  the  thickness) 

En  =  1 3.45x1 09  Pa  (through  thickness  modulus) 

E22  =  69.0x1 09  Pa  (inplane  modulus) 

v)2  =  0.04  (transverse)  and  v13  =  0.3  (inplane) 
damage  parameters  (see  section  2. 1  and  [22]) 

rij  =  l.OxlO-*,  n,  =  1.0 ,  oO0  =  70xl06  Pa  (3.5b) 

and  the  dispersion  parameters  (see  [22]) 

X,  =0,  X2=  5.0x10^,  y,  =  0,  y2=0.81  (3.5c) 

In  Figure  3.1,  a  WONDY  model  for  the  composite  plate  and  fluid  media  is  indicated.  For  these 
ID  analyses,  a  very  fine  finite  difference  grid  (zones)  was  employed  to  capture  the  damage  and  to 
accurately  model  the  shock  wave  as  it  moved  throught  the  fluid  and  into  the  plate.  For  full  2D 
(and  especially  3D  problems),  which  are  discussed  in  section  3.2,  a  much  coarser  grid  must 
obviously  be  employed,  otherwise  the  models  would  become  far  too  large  to  efficiently  analyze. 
Also,  we  note  in  Figure  3. 1  that  only  1  inch  of  the  fluid  was  included  in  the  model,  and  no  silent 
or  non-reflecting  boundary  conditions  were  required  at  point  A.  In  this  model,  our  concern  was 
strictly  with  the  damage  that  would  occur  once  the  compressive  wave  moves  through  the  fluid 
into  the  plate  and  is  reflected  at  the  free  end  B  as  a  tensile  wave.  The  analysis  was  terminated 
before  any  spurious  reflections  from  A  could  reach  the  plate. 
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Figures  3.2a,b  and  c  summarize  the  WONDY  damage  predictions  for  2,  4  and  8  inch  thick 
graphite/peek  plates.  The  estimated  boundaries  of  incipient  and  complete  spallation  are  indicated 
in  these  three  plots  in  which  range  and  weight  of  charge  have  been  varied.  We  note  that  in  these 
WONDY  analyses  dispersion/viscoelastic  effects  have  been  included  along  with  the  ID  damage 
models  (see  [22]).  Including  dispersion/viscoelastic  material  behavior  tends  to  slightly  mitigate 
the  overall  damage  predictions  due  to  the  material  damping  that  is  introduced  in  the  structural 
model. 

Next  we  will  apply  the  2D  transversely  isotropic  continuum  damage  theory  [19-20],  which  has 
been  programmed  into  ABAQUS/EXPLICIT  to  the  ID  fluid-structure  problem  indicated  in 
Figure  3.1.  (See  section  2.1  for  a  discussion  of  the  2D  damage  theory,  and  [22]  for  the  remaining 
values  of  the  various  elastic  variables  and  damage  parameters  used  for  the  graphite/peek  plate). 
The  2D  damage  theory  does  not  explicitly  contain  viscoelastic  capabilities  that  can  account  for  all 
the  damping  that  is  present  in  these  types  of  composite  structures.  Figure  3.3  describes  the  2D 
model  for  this  ID  problem.  Note  that  for  this  2D  model,  a  pseudo  infinite  element  and  several 
buffer  fluid  elements  were  employed  to  prevent  any  reflections  from  the  top  of  the  fluid  model.  In 
general,  this  was  not  necessary  since  the  analysis  was  terminated  before  any  reflections  from  the 
top  of  the  fluid  could  reach  the  plate.  A  total  of  156  2D  elements  were  used  -  55  in  the  fluid  and 
100  through  the  plate.  The  FEM  model  is  thus  relatively  coarse  compared  to  the  very  fine 
WONDY  model  indicated  in  Figure  3.1,  but  it  will  still  be  sufficient.  Figure  3.4  indicates  the 
through  thickness  damage  profile  predictions  from  the  ABAQUS/EXPLICIT  model  for  a  2  inch 
thick  graphite/peek  plate  subjected  to  a  charge  of  10  kg  at  a  range  of  one  meter.  In  general,  some 
damage  is  indicated  near  the  top  of  the  plate  in  Figure  3  .4,  but  there  is  no  spallation.  This  results 
compares  reasonably  well  to  the  through  thickness  damage  results  from  a  refined  WONDY 
analysis  (not  shown). 

Overall,  the  WONDY  and  ABAQUS/EXPLICIT  models  that  were  run  indicate  that  based  on 
ID  models  the  charges  need  to  be  very  close  in  to  do  significant  delamination  damage  to  the 
composite  graphite/peek  plates  considered. 
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Figure  3.1  WONDY  model  for  a  graphite/peek  plate  subjected  to  underwater 
shock  (ID  analysis). 
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Weight  of  charge 


Figure  3.2a  Predicted  spallation  (VI)  damage  for  TNT  in  water  and  a  graphite/ 
peek  plate  (2”). 

Weight  of  charge 
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Weight  of  charge 


Figure  3.2c  Predicted  spallation  (VI)  damage  for  TNT  in  water  and  a  graphite/ 
peek  plate  (8”). 
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Figure  3.3  ABAQUS/EXPLICIT  2D  model  for  ID  fluid-structure 
interaction. 
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Figure  3.4  Damage  results  through  a  2"  plate  from  an  ABAQUS/EXPLICIT  fluid- 
structure  interaction  model. 


4.  3D  CONTINUUM  DAMAGE 

The  development  of  a  suitable  3D  continuum  damage  theory  for  transient  problems  has  not  been 
an  easy  undertaking.  A  simple  extension  to  the  2D  theory  [19-21]  was  not  viable.  The 
heterogeneity  of  composite  structures,  especially  in  the  through  thickness  direction,  as  well  as  the 
directionality  of  the  damage  requires  some  sort  of  additional  homogenization  assumptions  to 
reduce  the  problem  to  a  managable  state. 

In  [22],  the  Nemes  and  Randles  [19-21]  2D  continuum  damage  model  (CDM)  was  extended  to 
3D.  This  section  is  essentially  a  brief  summary  of  the  3D  phenomenological  theory  or  CDM3D 
that  is  contained  in  [22].  At  the  time  of  the  writing  of  this  report,  a  preliminary  version  of  the 
CDM3D  has  been  programmed  into  ABAQUS/EXPLICIT  as  a  user  written  subroutine,  but  this 
capability  has  not  yet  been  tested. 

Essentially,  the  2D  continuum  damage  model  assumes  transverse  isotropy  in  the  1-2  plane  as 
indicated  in  Figure  4.1.  (Note  the  relabeling  of  the  1,2,  and  3  axes  in  Figure  4. 1  from  Figure  2.1.1 
to  allow  the  3  direction  to  be  through  the  thickness  of  the  composite).  Delamination  damage, 
representing  a  network  of  matrix  cracks  aligned  in  the  1-2  plane  with  normals  in  the  the  3 
direction  is  denoted  by  V3  (V,  in  [19-21])  and  depicted  in  Figure  4.2.  Inplane  matrix  cracks, 

which  may  traverse  or  lie  between  reinforcing  fibers  but  not  break  fibers,  are  indicated  in  Figure 
4.3  and  referred  to  as  Vs  damage  in  [19-21],  In  this  2D  approach,  damage  effects  were  tracked 

and  applied  to  the  engineering  moduli,  shear  moduli,  and  Poisson’s  ratios. 

The  assumption  of  transverse  isotropy  in  the  2D  formulation  represents  a  reasonable 
approximation  that  allows  the  through  thickness  direction  to  be  more  easily  discretized.  This 
approximation  eliminates  the  requirement  of  tracking  each  individual  ply  and  the  various 
compliances  associated  with  them.  It  represents  essentially  a  homogenization  of  the  laminate. 
Note  ,  however,  a  large  number  of  elements  are  still  required  in  the  through  thickness  direction  to 
accurately  model  transient  through  thickness  effects. 

In  a  3D  model  of  a  laminated  composite,  the  transverse  isotropy  assumption  is  no  longer  valid. 
Total  damage  which  will  be  represented  by  the  vector  D ,  can  be  expressed  is  the  same  form  as 
(2.1)  or: 

D  =  D,e, +D2e2+D3e3  (4.1) 
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where  again  e ;  are  the  unit  base  vectors  and  D,  are  individual  components  of  damage  with 
normals  aligned  in  the  e,  direction. 

A  direct  extension  of  the  2D  damage  model  to  3D  seems  logical  at  first  glance.  (In  the  2D 
theory  D,  and  D2or  V2  and  V3  in  section  2  and  [19-21]  are  simply  combined  into  Vs). 

However,  a  direct  extension  from  2D  is  very  inefficient  because  it  would  necessitate  the  tracking 
of  each  individual  ply  and  its  compliances. 

An  alternate  approach  to  the  modeling  of  3D  damage  is  to  track  the  damage  and  apply  it  to  the 
stresses  as  is  done  in  [32-33],  as  opposed  to  the  compliances  as  is  done  in  the  2D  model.  Solving 
(1.1)  for  D  yields: 

D  =  ^  (4.2) 

A 

D  can  thus  be  interpreted  as  the  relative  reduction  in  area  caused  by  damage  due  to  micro- 
cracking.  For  general  anisotropic  damage  in  three  dimensions,  the  effective  stress  a  can  be 
expressed  in  a  tensorial  form  as: 


a  =  M(D)o 


(4.3) 


where  Mis  the  material  damage  matrix,  and  Dis  the  damage  vector  whose  components  are 
defined  in  a  fashion  similar  to  (4.2).  In  [32-33],  the  following  form  of  M  is  proposed  to  account 
for  anisotropic  material  damage  in  the  principal  coordinate  system: 


M  = 
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(4.4) 


where,  as  in  2D  damage  theory,  the  fractions  (Ka^l,  I=l,2,3  have  been  added  to  prevent  a 
complete  loss  of  material  integrity  as  the  individual  components  of  damage  D,  go  to  1. 

For  individual  plies  (4.3)  and  (4.4)  can  be  easily  applied,  since  the  principal  directions  are  in  the 
1-2  plane  (see  Figure  4.1)  and  perpendicular  to  it.  Then  the  stress  tensor  a  can  be  transformed  to 
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the  global  1-2  axes.  In  this  work,  we  will  apply  (4.3)  and  (4.4)  to  groups  of  plies,  and  assume  the 
additional  approximation  that  the  global  1-2  axes  represent  the  principal  axes.  In  general  for 
laminated  composites,  the  1  and  2  directions  will  not  correspond  to  the  principal  directions.  This 
assumption,  however,  can  be  partially  rectified  by  adjusting  the  thresholds  for  damage  to  more 
accurately  reflect  the  grouping  of  the  individual  plies  within  a  single  finite  element. 

For  linear  elasticity,  the  relationship  between  the  stress  and  the  strain  tensors  can  be  written  as: 

5  =  Ce  (4.5) 

where  C  is  the  constitutive  tensor  and  8  is  the  strain.  In  the  case  of  an  orthotropic  lamina  in 
which  there  is  no  coupling  between  the  shear  stresses: 
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where  1,2  and  3  are  the  principal  axes.  In  terms  of  the  engineering  constants,  the  components  of 
the  constitutive  tensor  C  are: 


and 


C„  =(l-v23v23)/(E2E3p)  (4.7) 

C12  —  (V12  +  v32v13)  /  (E,E3P)  =  c21 
C)3  —  (v13  +  v12v23)  /  (E,E2P)  =  C3) 

C„«0-v„v,1)/(E1EJW 

^23  ~  (V23  +  V21  I  (EjEjP)  =  ^32 

CM=(l-vl2vJ,)/(E,E,P) 

C44=G12,  ^55  =  Gjj  =  G23 


3  =  0-  v12v21  -  v23v32  -  v31v13  -  2v21v32v13)/  (E,E2E3) 
where  v# ,  E{  and  G,  represent  the  various  Poisson  ratios,  the  engineering  and  shear  moduli.  In 
addition,  there  is  also  the  following  relationship: 
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(48) 


v  IE-  =  v  /E 

1J  I  J1  J 

in  which  the  summation  convention  has  been  suspended. 

Equations  (4.3)  -  (4.7)  can  be  constructed  for  each  individual  ply.  However  as  discussed 
previously,  it  is  not  our  intent  to  track  the  damage  in  each  ply  but  rather  in  groups  of  plies.  This 
latter  approach  is  much  more  feasible.  For  3D  and  2D  solid  finite  element,  explicit  transient 
analysis,  programs  such  as  ABAQUS/EXPLICIT  employ  linear  (interpolation)  solid  elements 
with  reduced  or  one  point  integration.  With  this  in  mind,  equations  (4.5)  -  (4.8)  can  be  applied  to 
groups  of  plies  and  their  homogenized  effect  into  one  or  more  elements.  With  a  linear  element 
due  to  one  point  integration,  the  stress  and  strain  do  not  vary  in  the  through  thickess  direction  (or 
in  the  other  two  direction  as  well).  Thus  (4.6)  for  a  group  of  N  plies  can  be  written  as: 

S=ZC.e  (4  9) 

1=1 

In  lumping  together  groups  of  plies  and  orientations  (a,b,c,d)  for  instance,  we  appear  to  be 
restricting  the  number  of  elements  allowed  in  the  3  direction.  In  some  instances  over  100 
elements  through  the  thickess  may  be  required  to  track  the  wave  nature  of  the  response  and 
ensuing  damage.  If  there  are  only  10  groups  of  plies  with  orientations  of  say  (a,b,c,d)  with  this 
approach  only  40  elements  in  the  3  direction  could  be  allowed  before  individual  plies  would  have 
to  be  modeled  with  more  than  one  element.  This  approach  would  also  be  very  difficult  with 
respect  to  pre-processing.  To  overcome  this  defficiency,  an  additional  approximation  is 
introduced,  which  assumes  each  group  of  plies  can  be  broken  down  into  subgroups  with  exactly 
the  same  group  of  plies  as  the  original  group  but  with  scaled  thicknesses.  Thus,  if  as  in  this  case, 
there  are  10  groups  of  plies  with  orientations  of  (a,b,c,d)  and  100  elements  through  the  thickness 
are  required  ,  is  assumed  that  there  are  100  groups  of  plies  (a,b,c,d)  with  thickness  1/10  of  the 
original  layup.  This  latter  homogenization  approximation  appears  to  be  a  reasonable  compromise 
to  the  modeling  of  each  individual  ply,  or  the  possibility  of  requiring  more  than  one  element  per 
ply  should  additional  elements  in  the  3  direction  be  necessary.  Of  course,  a  more  drastic 
homogenization  simplification  such  as  assuming  averaged  orthotropic  properties  through  the 
entire  thickness  can  also  be  made  with  equations  (4.2)  -  (4.7). 

We  now  return  to  (4.3)  and  (4.4).  The  strain  energy  W  of  an  undamped  material  can  be 
expressed  as: 
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W  =  — aTe 
2 


(4.10) 


Using  (4.6),  (4.10)  becomes: 

W  =  j6tCs  (4.11) 

in  terms  of  the  strain  tensors  or 

W  =  Vc-'g  (4.12) 

with  respect  to  the  stress  tensor.  The  strain  energy  WD  for  the  damaged  material  is  of  the  form: 

WD=^iTC''a  (4.13) 

Introducing  (4.3)  produces: 

wd=^5t(MtC-1M)o 

or 

WD=^5TC'1g  (4.14) 

where 

C'1  =MTC*'M  (4.15) 

Inverting  (4.15)  produces  the  damaged  constitutive  tensor: 

C  =  M-,C=(Mt)-'  (4.16) 

Because  of  (4.4),  M  is  assumed  to  be  a  diagonal  matrix  and  thus  its  inverse  is  easy  to  obtain. 

Having  postulated  an  approach  to  model  3D  continuum  damage,  the  specific  details  pertaining 
to  the  evolving  3D  damage  model  are  described.  The  intended  application  is  for  composite 
materials  which  usu sally  contain  a  great  number  of  imperfection  sites  from  which  cracks  can 
grow.  Thus,  the  nucleation  of  new  cracks  is  not  addressed.  Rather,  the  focus  is  on  the  growth 
and  evolution  of  damage. 
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Recall  that  the  damage  D  is  assumed  to  be  a  vector  with  the  components  D,  in  the  1,2  and  3 
directions.  As  in  the  2D  model  [19-20],  the  evolution  of  individual  components  of  damage  are 
assumed  to  be  governed  by  a  threshold  in  the  form: 

F,  (g,  f,  (Dj ))  <  0  for  no  damage  growth  (4.17) 

>  0  for  damage  growth 

where  (ij)=  1,2,3  and  there  is  no  sum  on  i;  F;  are  scalar  threshold  functions;  g  is  the  current 
stress  tensor;  and  f;  is  an  array  of  current  threshold  parameters  which  are  functions  of  Dj ,  the 
damage  component.  F(  is  assumed  to  be  of  the  Mohr-Coulomb  yield  surface  [28]  type  and  to  be 
dependent  upon  the  scalar  stress  measures  a,  and  x  (tension  and  shear)  as  follows: 

F,(ot,T,f)  =  (l+<7L)i),'!  -(f„  -o,)/f„  (4,18) 

*3i 

In  (4.18),  the  parameters  are  related  to  specific  growth  threshold  strengths  o0iand  xGi,  and  the 
Coulomb  friction  tangent  cpGi  as: 

°Gi  =  fli  —  ^2i 

tGi=4((fll/f2l)2-l)1/2  (4.19) 

^Gi  =  fji  /  ^2i 

The  tension  and  shear  growth  thresholds  as  well  as  the  Coulomb  friction  tangent  in  (4.19)  are 
assumed  to  be  functions  of  the  damage  D, .  The  threshold  surface  (F=0)  defined  by  (4.18)  is 
shown  in  Figure  4.4  along  with  the  various  parameters.  Note  that  in  this  figure  d  is  the  shortest 
distance  from  an  external  state  (ot,  x)  to  the  threshold  surface. 

Using  (4.18)  and  (4.19),  we  further  develop  the  details  for  the  components  of  damage  D;  and 
its  evolution.  For  D3  type  of  damage,  which  corresponds  to  delamination,  the  stress  measures 
c, i  and  x,  are  postulated  to  be  of  the  form: 

ot3=o33  (4.20) 

*3  “  (^31  +®32  ) 

The  growth  threshold  stresses  and  the  Coulomb  friction  tangent  in  (4. 19)  are  then  taken  to  be: 

°G3  =  0  “  D32  )oG30 
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(4.21) 


^G3  —  0  ^3  )^G30 

^G3  —  ^G30  ^3  (*Pg31  ~  *Pg3o) 

where  the  “0”  subscript  denotes  the  initial  undamaged  properties.  Note  that  for  oG3  and  xG3  all 
resistance  to  damage  is  lost  as  D3  goes  to  1;  and  the  Coulomb  friction  tangent  <pG1  can  vary 
linearly  with  D32  as  the  damages  progresses. 

From  (4.18)  and  (4.19),  the  variables  f23,  f13  and  f33  are: 


1  T  ^ 
f  -  (  G3 

23  2  %G12oG3 


^13  ~  °G3  +  ^23 


-^G3) 


(4.22) 


^33  “  G3 

To  complete  the  continuum  damage  formulation  for  D3  damage,  evolution  equations  are 
required.  From  [19-21],  the  following  relationship  is  postulated  for  D3  (delamination)  damage: 


F3(d3,D3) 


(4.23) 


where  d3  is  the  shortest  distance  in  the  at3,  t3  stress  space  -  see  Figure  4.4  for  the  threshold 

surface  F3  =  0.  For  d3  =  0,  F3  =  0  and  — —  =  0  for  all  stress  points  interior  or  on  the  threshold 

dt 

surface.  The  specific  form  of  (4.23)  employed  for  the  3D  theory  is  similar  to  (2.9a)  with  D3 
replacing  V,  and  the  subscript  3  replacing  1  or: 


-jP  =  (d3  /  <*G30  )"3  /  (Tl3  (1  -  D32 ) 


(4.24) 


Next  we  consider  D,  and  D2  type  of  damage,  which  relate  to  matrix  cracking.  The  discussion 
will  focus  on  D, ,  but  D2  type  of  damage  will  be  assumed  to  be  of  the  same  form.  For  D,  type  of 
damage,  the  stress  measures  otl  and  t,  required  for  (4. 18)  and  (4. 19)  are  postulated  to  be  of  the 
form: 


^ti  “  On 


x,  =(o122+o132)1 


(4.25) 
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Equation  (4.25)  are  used  in  conjunction  with  (4  .18)  and  (4. 1 9).  The  form  of  the  growth  threshold 
stresses  is  not  of  the  same  form  as  the  D3  damage,  which  is  given  by  (4.20).  D,  (and  also  D2 ) 

damage  consist  of  an  ever  denser  network  of  cracks  which  tend  to  saturate.  There  is  not 
acceleration  to  a  castrophic  failure  as  in  the  case  of  D3  damage  which  represents  delamination. 

Saturation  can  be  forced  through  the  growth  thresholds  as  done  in  [19-20]  or: 

°g,=<W0-D,2)  (4.26) 

^G1  =  *G10  !  0  —  Dj  ) 

while  the  Coulomb  friction  tangent  is  taken  to  be 

^Pgi  =  ^Pgio (^Pgii  ^Pgio)  (4.27) 

In  (4.26),  it  is  seen  that  as  D,goes  to  1  damage  becomes  more  difficult  since  the  growth 
thresholds  aG1  and  tG)  greatly  increase. 

For  D,  (and  D2 )  damage,  the  final  ingredient  required  is  the  specific  form  of  the  rate  evolution 
equation  described  by  (4.23).  Since  there  is  no  acceleration  to  catastrophic  failure,  a  form  similar 
to  (2.9b)  for  V,  damage  in  the  2d  theory  is  employed: 

^jL  =  (d1/o<11. )-/!,,  (4.28) 

where  as  in  (2.9b),  t),  is  a  time  constant,  d,  is  normalized  with  the  virgin  growth  threshold  stress 
oG10 ,  and  n,  is  a  positive  term  for  the  dimensionless  stress  distance. 

Another  mode  of  continuum  damage  is  the  breakage  of  the  reinforcing  fibers.  To  account  for 
this  effect,  a  simple  maximum  strain  criteria  is  used  in  the  3D  model.  The  D, ,  D2  and  D3  forms 

of  damage  cause  softening  of  the  composite  response,  and  thus  directly  influence  this  maximum 
strain  criteria. 
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60  deg  ply  layup. 


FIGURE  4.4 


Threshold  surface  for  the  onset  of  damage. 
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5.  SUMMARY  AND  CONCLUSIONS 

In  this  report,  we  have  discussed  in  the  Introduction  some  of  the  difficulties  associated  with  the 
modeling  of  the  shock  response  and  damage  in  marine  composites.  The  existing  technical  barriers 
are  indeed  formidable  but  some  progress  has  been  made.  Several  key  areas  that  must  be 
addressed  are:  dispersion  and  viscoelasticity  in  the  composite  structure;  homogenization 
assumptions;  through  thickness  modeling  and  damage;  and  dynamic  versus  static  properties  and 
behavior. 

In  sections  2-5,  evolving  continuum  damage  models  based  primarily  on  the  work  of  Nemes 
and  Randles  [19-21]  have  been  discussed  with  application  to  ID  and  2D  impact  and  underwater 
shock  problems.  This  approach,  which  has  been  extended  into  3D  in  [22],  is  phenomenological  in 
nature  requiring  the  determining  of  several  empirical  constants  for  specific  composite  structures. 
Comparisons  to  experimental  data  in  [23]  are  good.  ID  parametric  studies  are  discussed  in 
sections  3.1  and  3.2  for  underwater  shock  loading  environments. 

In  general,  the  phenomenological  approach  taken  here  at  NRL  toward  the  modeling  of  the 
shock  response  and  damage  in  composites  is  of  utility  but  is  of  course  limited.  The  first  author 
likens  this  situation  to  turbulence  modeling  in  fluid  dynamics  [34-35],  No  general  purpose 
turbulence  theory  exists.  The  current  turbulence  models  are  applicable  to  specific  flow  situations 
and  usually  require  some  emphircal  constants  to  be  used  to  obtain  accurate  simulations.  The 
statistical  nature  of  turbulence  is  addressed  using  time  averaging  techniques  and  generating 
continuum  equations  containing  fluctuations. 

As  with  turbulence,  the  material  properties  as  well  as  the  shock  response  and  evolving  damage 
in  composite  structures  have  a  non-deterministic,  statistical  nature  that  is  difficult  to  adequately 
address  with  continuum  mechanics  analysis  methods.  The  continuum  assumption  eventually 
either  breaks  down  or  must  be  modified  as  one  goes  from  the  macro  to  meso  to  micro  scales.  The 
interaction  of  these  time  scales  is  difficult  to  incorporate  into  a  macro  finite  element  model  that  is 
robust  and  accurate  for  a  wide  variety  of  materials  and  structures. 
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